一.从FNA文件中提取指定区间的序列

方法1：使用seqkit工具（推荐）

# 安装seqkit (如果尚未安装)
conda install -c bioconda seqkit

# 提取指定区间 (例如：从100到200的位置)
seqkit subseq --chr "序列ID" -r 100:200 input.fna > output.fna

方法2：使用samtools faidx

# 首先建立索引
samtools faidx input.fna

# 然后提取区间
samtools faidx input.fna "序列ID:100-200" > output.fna

方法3：使用Python脚本

from Bio import SeqIO

def extract_sequence(fna_file, seq_id, start, end, output_file):
    for record in SeqIO.parse(fna_file, "fasta"):
        if record.id == seq_id:
            subseq = record.seq[start-1:end]  # Python是0-based，生物坐标是1-based
            with open(output_file, "w") as f:
                f.write(f">{record.id}_{start}_{end}\n{subseq}\n")
            return
    print(f"序列ID {seq_id} 未找到")

# 使用示例
extract_sequence("input.fna", "chr1", 100, 200, "output.fna")

方法4：使用bedtools

# 首先创建一个BED文件定义区间
echo -e "序列ID\t99\t200" > region.bed  # BED是0-based半开区间

# 然后提取序列
bedtools getfasta -fi input.fna -bed region.bed -fo output.fna

注意事项

1. 区间定义：不同工具对区间定义可能不同(0-based或1-based，闭区间或半开区间)
2. 序列ID：确保输入的序列ID与文件中的完全匹配
3. 大文件处理：对于大基因组文件，建议先建立索引(如samtools faidx)
4. 反向互补链：如果需要提取反向互补链，可以添加参数如--reverse-complement(seqkit)或处理Python中的序列

二.从FASTA文件 map_gene.fa 中保留最长的200条序列

方法2：使用 seqkit（推荐）

seqkit sort --by-length --reverse map_gene.fa \
  | seqkit head -n 200 > top200.fa

步骤解释：

1. seqkit sort 按序列长度降序排序。
2. seqkit head 提取前200条记录。

安装 seqkit：

wget https://github.com/shenwei356/seqkit/releases/download/v2.3.0/seqkit_linux_amd64.tar.gz
tar -zxvf seqkit_linux_amd64.tar.gz
sudo mv seqkit /usr/local/bin/

方法3：纯 awk + 基础命令

awk '/^>/ { 
        if (header) print header ORS seq
        header=$0
        seq=""
        next
     } 
     { seq = seq $0 } 
     END { if (header) print header ORS seq }' map_gene.fa \
  | awk 'BEGIN {RS=">"; FS="\n"} NR>1 { 
        len=0; seq=""
        for (i=2; i<=NF; i++) { len+=length($i); seq=seq $i }
        print len, ">"$1, seq
     }' \
  | sort -k1,1nr \
  | head -n 200 \
  | cut -d' ' -f2- \
  | tr ' ' '\n' \
  | awk '/^>/ {print; next} {gsub(/.{80}/,"&\n"); printf "%s", $0}' > top200.fa

步骤解释：

1. 合并多行序列为单行。
2. 计算每条序列的长度并附加到行首。
3. 按长度降序排序，取前200条。
4. 恢复FASTA格式，每80字符换行。

三.批量修改FASTA文件序列ID

方法1：使用Linux命令 awk

awk 'BEGIN {
    # 加载chrom_map.txt到关联数组map中
    while (getline < "chrom_map.txt") {
        map[$1] = $2
    }
}
# 处理以">"开头的行（序列ID行）
/^>/ {
    # 提取旧ID（去掉开头的">"，取第一个字段）
    old_id = substr($1, 2)
    # 如果旧ID在映射表中，则替换为新ID
    if (old_id in map) {
        $1 = ">" map[old_id]  # 修改第一个字段为新ID
    }
}
# 打印所有行（1表示默认动作：打印当前行）
1' input.fasta > output.fasta

四.从FASTA格式基因组中批量提取指定染色体

方法一：使用 samtools 工具（推荐）

适用场景：处理大文件时效率高，需安装 samtools。

1. 生成索引文件：

   samtools faidx genome.fasta
   

   这会生成 genome.fasta.fai 索引文件。

2. 提取指定染色体：

   samtools faidx genome.fasta chr1 chr3 chr5 > output.fasta
   

   若染色体列表存储在文件（如 chr_list.txt）中，可使用：

   xargs samtools faidx genome.fasta < chr_list.txt > output.fasta
   

方法二：使用 awk 命令

适用场景：快速提取，适合简单需求。

awk -v targets="chr1 chr3 chr5" '
BEGIN {
  split(targets, arr, " ");
  for (i in arr) chrs[arr[i]] = 1
}
/^>/ {
  chr_name = substr($1, 2);  # 提取">"后的名称
  if (chr_name in chrs) {
    print;
    flag = 1;
    next;
  }
  flag = 0;
}
flag { print }
' genome.fasta > output.fasta

方法三：使用 Python 脚本

from Bio import SeqIO

input_file = "genome.fasta"
output_file = "output.fasta"
target_chromosomes = {"chr1", "chr3", "chr5"}  # 修改为目标染色体名称

# 读取并筛选记录
records = []
for record in SeqIO.parse(input_file, "fasta"):
    if record.id in target_chromosomes:
        records.append(record)

# 写入输出文件
SeqIO.write(records, output_file, "fasta")